rm(list = ls())
library(data.table)
library(foreign)
library(kernlab)
library(MASS)  # for mvrnorm
library(ggplot2)
library(gridExtra)
library(glmnet)
library(caret)
library(memisc) #cases
setwd('~/github/bdr/')
source('functions.R')

Read in data

# import data
data_sept18 = data.table(read.spss('data/Sept18/Sept18 public.sav', to.data.frame = T), stringsAsFactors = F)
Undeclared level(s) 2, 3, 4 added in variable: densityUndeclared level(s) 2, 3, 4 added in variable: sdensityUndeclared level(s) 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95 added in variable: ageUndeclared level(s) 1, 2, 3, 4, 5, 6, 7 added in variable: hh1Undeclared level(s) 1, 2, 3, 4, 5, 6, 7 added in variable: hh3
data_sept18
data_june18 = data.table(read.spss('data/June18/June18 public.sav', to.data.frame = T), stringsAsFactors = F)
Undeclared level(s) 2, 3, 4 added in variable: densityUndeclared level(s) 2, 3, 4 added in variable: sdensityUndeclared level(s) 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94 added in variable: ageUndeclared level(s) 1, 2, 3, 4, 5, 6, 7 added in variable: hh1Undeclared level(s) 1, 2, 3, 4, 5, 6 added in variable: hh3
data_june18
data_may18 = data.table(read.spss('data/May18/May18 public.sav', to.data.frame = T), stringsAsFactors = F)
Undeclared level(s) 2, 3, 4 added in variable: densityUndeclared level(s) 2, 3, 4 added in variable: sdensityUndeclared level(s) 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 15, 20, 24 added in variable: slepnhisUndeclared level(s) 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 92, 94, 95 added in variable: ageUndeclared level(s) 1, 2, 3, 4, 5, 6, 7 added in variable: hh1Undeclared level(s) 1, 2, 3, 4, 5, 6, 7 added in variable: hh3
data_may18
# data_march18 = data.table(read.spss('data/March18/March18 public.sav', to.data.frame = T), stringsAsFactors = F)
# data_march18
###### PREP DATA
#specify which should be survey
survey_vars = c('demo_mode', 'demo_education', 'demo_phonetype', 'month_called', 'demo_ideology')
file_and_survey_vars = c('demo_sex', 'demo_age_bucket', 'demo_state', 'demo_income', 'demo_region', 'demo_race', 'demo_hispanic')
# covars = c('sample', 'int_date', 'sex', 'educ', 'hisp', 'racecmb', 'relig', 'income', 'party', 'hh1', 'ql1', 'sstate', 'sdensity')
# 
# # check that covariates are in all data sets
# covars %in% names(data_june18)
# covars %in% names(data_sept18)
# covars %in% names(data_may18)
# covars %in% names(data_march18)
data_sept18
data_may18[, .N, party]

Recode Data

NAs introduced by coercionconditions are not mutually exclusivecondition is.na(age_num) is never satisfiedconditions are not mutually exclusiveconditions are not mutually exclusiveconditions are not mutually exclusiveconditions are not mutually exclusiveconditions are not mutually exclusiveconditions are not mutually exclusiveconditions are not mutually exclusiveconditions are not mutually exclusiveAdding new column 'qsupport' then assigning NULL (deleting it).NAs introduced by coercionconditions are not mutually exclusivecondition is.na(age_num) is never satisfiedconditions are not mutually exclusiveconditions are not mutually exclusiveconditions are not mutually exclusiveconditions are not mutually exclusiveconditions are not mutually exclusiveconditions are not mutually exclusiveconditions are not mutually exclusiveconditions are not mutually exclusiveAdding new column 'qsupport' then assigning NULL (deleting it).NAs introduced by coercionconditions are not mutually exclusivecondition is.na(age_num) is never satisfiedconditions are not mutually exclusiveconditions are not mutually exclusiveconditions are not mutually exclusiveconditions are not mutually exclusiveconditions are not mutually exclusiveconditions are not mutually exclusiveconditions are not mutually exclusiveconditions are not mutually exclusiveAdding new column 'qsupport' then assigning NULL (deleting it).

Categorize variables

# create data table with vars and levels
covars = names(data_recoded)[grepl('demo', names(data_recoded))]
covars = data.table(do.call(rbind, lapply(covars, function(c){
  cbind(c, data_recoded[, .(level = unique(get(c)))][order(level)])
})))
setnames(covars, c('var', 'level'))
covars[, level_modmat := paste0(var, level)]
covars[, in_survey := as.numeric(var %in% survey_vars)]
covars[, in_both := as.numeric(var %in% file_and_survey_vars)]
covars[, in_file := as.numeric(in_survey + in_both == 0)]

Get data partitions

Proability of surveyed

# scale age and set NAs to 0
data_recoded[, age_scaled := scale(age_num)]
data_recoded[is.na(age_scaled), age_scaled := 0]
data_recoded[, p_surveyed := 
               (-2)
             + age_scaled 
             - 0.5 * is.na(age_num) 
             + 1.5 * as.numeric(demo_mode == 'cell') 
             - 1.5 * as.numeric(demo_party == '05-Ind')
             - 3 * as.numeric(demo_party == "99-DK/refused") 
             + 1.5 * as.numeric(demo_education %in% c('01-postgrad', '02-bach'))
             + 3 * as.numeric(demo_ideology == 'Very conservative' | demo_ideology == 'Very liberal')
             ]
data_recoded[, p_surveyed := exp(p_surveyed)/(1 + exp(p_surveyed))]
hist(data_recoded[, p_surveyed])

data_recoded[, .(.N, mean(p_surveyed)), .(demo_age_bucket)][order(demo_age_bucket)]
data_recoded[, .(.N, mean(p_surveyed)), .(demo_mode)][order(demo_mode)]
data_recoded[, .(.N, mean(p_surveyed)), .(demo_party)][order(demo_party)]
data_recoded[, .(.N, mean(p_surveyed)), .(demo_ideology)][order(demo_ideology)]

Probability of matched

data_recoded[, p_matched := NULL]
Adding new column 'p_matched' then assigning NULL (deleting it).
data_recoded[, p_matched :=
               -2 +
               -2 * as.numeric(demo_mode == 'landline') 
             + 3 * as.numeric(demo_race == 'W')
             + -2 * as.numeric(demo_reg == '03-No')
             + -1 * as.numeric(demo_hhsize == 2)
             + -2 * as.numeric(demo_hhsize == 3)
             + 0.5*age_scaled
             + as.numeric(demo_income)/3
             - 4* as.numeric(demo_income == '99-DK/refused')
             ]
data_recoded[, p_matched := exp(p_matched)/(1 + exp(p_matched))]
hist(data_recoded$p_matched)

data_recoded[, .(.N, mean(p_matched)), demo_mode]
data_recoded[, .(.N, mean(p_matched)), demo_hispanic]
data_recoded[, .(.N, mean(p_matched)), demo_age_bucket][order(demo_age_bucket)]

Test/training sets

testtrain = getTestTrain(data = data_recoded
             , n_holdout = 1000, n_surveyed = 2000, n_matched = 1000
             , p_surveyed = data_recoded$p_surveyed
             , p_matched = data_recoded$p_matched
             )
Adding new column 'holdout' then assigning NULL (deleting it).Adding new column 'surveyed' then assigning NULL (deleting it).Adding new column 'matched' then assigning NULL (deleting it).
data_recoded[, .N, list(holdout, surveyed, matched, voterfile)]
data_recoded = testtrain$data
data_recoded[, .(.N, prop_surveyed = mean(surveyed)), demo_age_bucket][order(demo_age_bucket)]
data_recoded[, .(.N, prop_surveyed = mean(surveyed)), demo_party][order(demo_party)]
data_recoded[, .(.N, mean(y_dem)), .(holdout, surveyed, matched)]
data_recoded[, mean(y_dem)]
[1] 0.4864043

Simple Dist Regression

Assign bags

Get landmark points

landmarks = getLandmarks(data = data_recoded
             , vars = unique(covars[in_both == 1 | in_file == 1,]$var)
             , n_landmarks = 12
             , subset_ind = (data_recoded$voterfile == 1))
X_file = landmarks$X[data_recoded$voterfile == 1, ]
n_landmarks

[1] 12

Embed file in feature space

Choose scale parameter \(\sigma\) for RBF kernel - use median heuristic for now

sigma

[1] 1.30073e-05

Do basic Dist Reg

Hyperparameter tuning

names(fit_basicDR)

[1] "data"      "fit"       "landmarks" "bags"      "y_hat"     "mse_test" 
param_nlandmarks

 [1]   12   20   33   55   90  148  245  403  665 1097

LS0tCnRpdGxlOiAiQkRSIFBldyBEYXRhIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpgYGB7cn0Kcm0obGlzdCA9IGxzKCkpCmxpYnJhcnkoZGF0YS50YWJsZSkKbGlicmFyeShmb3JlaWduKQpsaWJyYXJ5KGtlcm5sYWIpCmxpYnJhcnkoTUFTUykgICMgZm9yIG12cm5vcm0KbGlicmFyeShnZ3Bsb3QyKQpsaWJyYXJ5KGdyaWRFeHRyYSkKbGlicmFyeShnbG1uZXQpCmxpYnJhcnkoY2FyZXQpCmxpYnJhcnkobWVtaXNjKSAjY2FzZXMKc2V0d2QoJ34vZ2l0aHViL2Jkci8nKQoKc291cmNlKCdmdW5jdGlvbnMuUicpCmBgYAoKCgojIyBSZWFkIGluIGRhdGEKYGBge3J9CgojIGltcG9ydCBkYXRhCmRhdGFfc2VwdDE4ID0gZGF0YS50YWJsZShyZWFkLnNwc3MoJ2RhdGEvU2VwdDE4L1NlcHQxOCBwdWJsaWMuc2F2JywgdG8uZGF0YS5mcmFtZSA9IFQpLCBzdHJpbmdzQXNGYWN0b3JzID0gRikKZGF0YV9zZXB0MTgKCmRhdGFfanVuZTE4ID0gZGF0YS50YWJsZShyZWFkLnNwc3MoJ2RhdGEvSnVuZTE4L0p1bmUxOCBwdWJsaWMuc2F2JywgdG8uZGF0YS5mcmFtZSA9IFQpLCBzdHJpbmdzQXNGYWN0b3JzID0gRikKZGF0YV9qdW5lMTgKCmRhdGFfbWF5MTggPSBkYXRhLnRhYmxlKHJlYWQuc3BzcygnZGF0YS9NYXkxOC9NYXkxOCBwdWJsaWMuc2F2JywgdG8uZGF0YS5mcmFtZSA9IFQpLCBzdHJpbmdzQXNGYWN0b3JzID0gRikKZGF0YV9tYXkxOAoKIyBkYXRhX21hcmNoMTggPSBkYXRhLnRhYmxlKHJlYWQuc3BzcygnZGF0YS9NYXJjaDE4L01hcmNoMTggcHVibGljLnNhdicsIHRvLmRhdGEuZnJhbWUgPSBUKSwgc3RyaW5nc0FzRmFjdG9ycyA9IEYpCiMgZGF0YV9tYXJjaDE4CgojIyMjIyMgUFJFUCBEQVRBCgojc3BlY2lmeSB3aGljaCBzaG91bGQgYmUgc3VydmV5CnN1cnZleV92YXJzID0gYygnZGVtb19tb2RlJywgJ2RlbW9fZWR1Y2F0aW9uJywgJ2RlbW9fcGhvbmV0eXBlJywgJ21vbnRoX2NhbGxlZCcsICdkZW1vX2lkZW9sb2d5JykKZmlsZV9hbmRfc3VydmV5X3ZhcnMgPSBjKCdkZW1vX3NleCcsICdkZW1vX2FnZV9idWNrZXQnLCAnZGVtb19zdGF0ZScsICdkZW1vX2luY29tZScsICdkZW1vX3JlZ2lvbicsICdkZW1vX3JhY2UnLCAnZGVtb19oaXNwYW5pYycpCgoKIyBjb3ZhcnMgPSBjKCdzYW1wbGUnLCAnaW50X2RhdGUnLCAnc2V4JywgJ2VkdWMnLCAnaGlzcCcsICdyYWNlY21iJywgJ3JlbGlnJywgJ2luY29tZScsICdwYXJ0eScsICdoaDEnLCAncWwxJywgJ3NzdGF0ZScsICdzZGVuc2l0eScpCiMgCiMgIyBjaGVjayB0aGF0IGNvdmFyaWF0ZXMgYXJlIGluIGFsbCBkYXRhIHNldHMKIyBjb3ZhcnMgJWluJSBuYW1lcyhkYXRhX2p1bmUxOCkKIyBjb3ZhcnMgJWluJSBuYW1lcyhkYXRhX3NlcHQxOCkKIyBjb3ZhcnMgJWluJSBuYW1lcyhkYXRhX21heTE4KQojIGNvdmFycyAlaW4lIG5hbWVzKGRhdGFfbWFyY2gxOCkKCmRhdGFfc2VwdDE4CgpkYXRhX21heTE4WywgLk4sIHBhcnR5XQpgYGAKCgojIyBSZWNvZGUgRGF0YQpgYGB7ciwgZWNobyA9IEYsIHdhcm5pbmc9Rn0KIyBSRUNPREUgYWxsIGRhdGEgc2V0cwpkYXRhX3JlY29kZWQgPSByYmluZGxpc3QobGFwcGx5KGxpc3QoZGF0YV9zZXB0MTgsIGRhdGFfbWF5MTgsIGRhdGFfanVuZTE4KSwgZG9QZXdSZWNvZGUpKQoKZGF0YV9yZWNvZGVkCmBgYAoKCiMjIENhdGVnb3JpemUgdmFyaWFibGVzCmBgYHtyfQojIGNyZWF0ZSBkYXRhIHRhYmxlIHdpdGggdmFycyBhbmQgbGV2ZWxzCmNvdmFycyA9IG5hbWVzKGRhdGFfcmVjb2RlZClbZ3JlcGwoJ2RlbW8nLCBuYW1lcyhkYXRhX3JlY29kZWQpKV0KY292YXJzID0gZGF0YS50YWJsZShkby5jYWxsKHJiaW5kLCBsYXBwbHkoY292YXJzLCBmdW5jdGlvbihjKXsKICBjYmluZChjLCBkYXRhX3JlY29kZWRbLCAuKGxldmVsID0gdW5pcXVlKGdldChjKSkpXVtvcmRlcihsZXZlbCldKQp9KSkpCnNldG5hbWVzKGNvdmFycywgYygndmFyJywgJ2xldmVsJykpCmNvdmFyc1ssIGxldmVsX21vZG1hdCA6PSBwYXN0ZTAodmFyLCBsZXZlbCldCgpjb3ZhcnNbLCBpbl9zdXJ2ZXkgOj0gYXMubnVtZXJpYyh2YXIgJWluJSBzdXJ2ZXlfdmFycyldCmNvdmFyc1ssIGluX2JvdGggOj0gYXMubnVtZXJpYyh2YXIgJWluJSBmaWxlX2FuZF9zdXJ2ZXlfdmFycyldCmNvdmFyc1ssIGluX2ZpbGUgOj0gYXMubnVtZXJpYyhpbl9zdXJ2ZXkgKyBpbl9ib3RoID09IDApXQpgYGAKCiMjIEdldCBkYXRhIHBhcnRpdGlvbnMKCiMjIyBQcm9hYmlsaXR5IG9mIHN1cnZleWVkCmBgYHtyfQojIHNjYWxlIGFnZSBhbmQgc2V0IE5BcyB0byAwCmRhdGFfcmVjb2RlZFssIGFnZV9zY2FsZWQgOj0gc2NhbGUoYWdlX251bSldCmRhdGFfcmVjb2RlZFtpcy5uYShhZ2Vfc2NhbGVkKSwgYWdlX3NjYWxlZCA6PSAwXQoKZGF0YV9yZWNvZGVkWywgcF9zdXJ2ZXllZCA6PSAKICAgICAgICAgICAgICAgKC0yKQogICAgICAgICAgICAgKyBhZ2Vfc2NhbGVkIAogICAgICAgICAgICAgLSAwLjUgKiBpcy5uYShhZ2VfbnVtKSAKICAgICAgICAgICAgICsgMS41ICogYXMubnVtZXJpYyhkZW1vX21vZGUgPT0gJ2NlbGwnKSAKICAgICAgICAgICAgIC0gMS41ICogYXMubnVtZXJpYyhkZW1vX3BhcnR5ID09ICcwNS1JbmQnKQogICAgICAgICAgICAgLSAzICogYXMubnVtZXJpYyhkZW1vX3BhcnR5ID09ICI5OS1ESy9yZWZ1c2VkIikgCiAgICAgICAgICAgICArIDEuNSAqIGFzLm51bWVyaWMoZGVtb19lZHVjYXRpb24gJWluJSBjKCcwMS1wb3N0Z3JhZCcsICcwMi1iYWNoJykpCiAgICAgICAgICAgICArIDMgKiBhcy5udW1lcmljKGRlbW9faWRlb2xvZ3kgPT0gJ1ZlcnkgY29uc2VydmF0aXZlJyB8IGRlbW9faWRlb2xvZ3kgPT0gJ1ZlcnkgbGliZXJhbCcpCiAgICAgICAgICAgICBdCmRhdGFfcmVjb2RlZFssIHBfc3VydmV5ZWQgOj0gZXhwKHBfc3VydmV5ZWQpLygxICsgZXhwKHBfc3VydmV5ZWQpKV0KaGlzdChkYXRhX3JlY29kZWRbLCBwX3N1cnZleWVkXSkKCmRhdGFfcmVjb2RlZFssIC4oLk4sIG1lYW4ocF9zdXJ2ZXllZCkpLCAuKGRlbW9fYWdlX2J1Y2tldCldW29yZGVyKGRlbW9fYWdlX2J1Y2tldCldCmRhdGFfcmVjb2RlZFssIC4oLk4sIG1lYW4ocF9zdXJ2ZXllZCkpLCAuKGRlbW9fbW9kZSldW29yZGVyKGRlbW9fbW9kZSldCmRhdGFfcmVjb2RlZFssIC4oLk4sIG1lYW4ocF9zdXJ2ZXllZCkpLCAuKGRlbW9fcGFydHkpXVtvcmRlcihkZW1vX3BhcnR5KV0KZGF0YV9yZWNvZGVkWywgLiguTiwgbWVhbihwX3N1cnZleWVkKSksIC4oZGVtb19pZGVvbG9neSldW29yZGVyKGRlbW9faWRlb2xvZ3kpXQpgYGAKCiMjIyBQcm9iYWJpbGl0eSBvZiBtYXRjaGVkCmBgYHtyfQpkYXRhX3JlY29kZWRbLCBwX21hdGNoZWQgOj0gTlVMTF0KZGF0YV9yZWNvZGVkWywgcF9tYXRjaGVkIDo9CiAgICAgICAgICAgICAgIC0yICsKICAgICAgICAgICAgICAgLTIgKiBhcy5udW1lcmljKGRlbW9fbW9kZSA9PSAnbGFuZGxpbmUnKSAKICAgICAgICAgICAgICsgMyAqIGFzLm51bWVyaWMoZGVtb19yYWNlID09ICdXJykKICAgICAgICAgICAgICsgLTIgKiBhcy5udW1lcmljKGRlbW9fcmVnID09ICcwMy1ObycpCiAgICAgICAgICAgICArIC0xICogYXMubnVtZXJpYyhkZW1vX2hoc2l6ZSA9PSAyKQogICAgICAgICAgICAgKyAtMiAqIGFzLm51bWVyaWMoZGVtb19oaHNpemUgPT0gMykKICAgICAgICAgICAgICsgMC41KmFnZV9zY2FsZWQKICAgICAgICAgICAgICsgYXMubnVtZXJpYyhkZW1vX2luY29tZSkvMwogICAgICAgICAgICAgLSA0KiBhcy5udW1lcmljKGRlbW9faW5jb21lID09ICc5OS1ESy9yZWZ1c2VkJykKICAgICAgICAgICAgIF0KZGF0YV9yZWNvZGVkWywgcF9tYXRjaGVkIDo9IGV4cChwX21hdGNoZWQpLygxICsgZXhwKHBfbWF0Y2hlZCkpXQpoaXN0KGRhdGFfcmVjb2RlZCRwX21hdGNoZWQpCgpkYXRhX3JlY29kZWRbLCAuKC5OLCBtZWFuKHBfbWF0Y2hlZCkpLCBkZW1vX21vZGVdCmRhdGFfcmVjb2RlZFssIC4oLk4sIG1lYW4ocF9tYXRjaGVkKSksIGRlbW9faGlzcGFuaWNdCmRhdGFfcmVjb2RlZFssIC4oLk4sIG1lYW4ocF9tYXRjaGVkKSksIGRlbW9fYWdlX2J1Y2tldF1bb3JkZXIoZGVtb19hZ2VfYnVja2V0KV0KYGBgCgojIyMgVGVzdC90cmFpbmluZyBzZXRzCmBgYHtyfQp0ZXN0dHJhaW4gPSBnZXRUZXN0VHJhaW4oZGF0YSA9IGRhdGFfcmVjb2RlZAogICAgICAgICAgICAgLCBuX2hvbGRvdXQgPSAxMDAwLCBuX3N1cnZleWVkID0gMjAwMCwgbl9tYXRjaGVkID0gMTAwMAogICAgICAgICAgICAgLCBwX3N1cnZleWVkID0gZGF0YV9yZWNvZGVkJHBfc3VydmV5ZWQKICAgICAgICAgICAgICwgcF9tYXRjaGVkID0gZGF0YV9yZWNvZGVkJHBfbWF0Y2hlZAogICAgICAgICAgICAgKQpkYXRhX3JlY29kZWRbLCAuTiwgbGlzdChob2xkb3V0LCBzdXJ2ZXllZCwgbWF0Y2hlZCwgdm90ZXJmaWxlKV0KCmRhdGFfcmVjb2RlZCA9IHRlc3R0cmFpbiRkYXRhCgpkYXRhX3JlY29kZWRbLCAuKC5OLCBwcm9wX3N1cnZleWVkID0gbWVhbihzdXJ2ZXllZCkpLCBkZW1vX2FnZV9idWNrZXRdW29yZGVyKGRlbW9fYWdlX2J1Y2tldCldCmRhdGFfcmVjb2RlZFssIC4oLk4sIHByb3Bfc3VydmV5ZWQgPSBtZWFuKHN1cnZleWVkKSksIGRlbW9fcGFydHldW29yZGVyKGRlbW9fcGFydHkpXQoKZGF0YV9yZWNvZGVkWywgLiguTiwgbWVhbih5X2RlbSkpLCAuKGhvbGRvdXQsIHN1cnZleWVkLCBtYXRjaGVkKV0KZGF0YV9yZWNvZGVkWywgbWVhbih5X2RlbSldCmBgYAoKIyMgU2ltcGxlIERpc3QgUmVncmVzc2lvbiAKCiMjIyMgQXNzaWduIGJhZ3MKYGBge3J9CiMjIyMgQ3JlYXRlIGJhZ3Mgb2Ygc3VydmV5IHJlc3BvbnNlcyB1c2luZyB2YXJpYWJsZXMgb2JzZXJ2ZWQgaW4gQk9USCBzdXJ2ZXkgQU5EIHZvdGVyZmlsZQojIHNob3VsZCB3ZSBtYWtlIHRoZSBiYWdzIHdpdGggdGhlIHN1YnNldCBvZiBkYXRhIGZyb20gdGhlIHN1cnZleSwgb3IgdGhlIHdob2xlIGZpbGU/Pz8KCmJhZ3MgPSBnZXRCYWdzKGRhdGEgPSBkYXRhX3JlY29kZWRbc3VydmV5ZWQgPT0gMSxdCiAgICAgICAgLCB2YXJzID0gZmlsZV9hbmRfc3VydmV5X3ZhcnMKICAgICAgICAsIG5fYmFncyA9IDMwCiAgICAgICAgLCBuZXdkYXRhID0gZGF0YV9yZWNvZGVkW3ZvdGVyZmlsZSA9PSAxIHwgaG9sZG91dCA9PSAxLCBdKQoKZGF0YV9yZWNvZGVkW3N1cnZleWVkID09IDEsIGJhZyA6PSBiYWdzJGJhZ3NdCmRhdGFfcmVjb2RlZFt2b3RlcmZpbGUgPT0gMSB8IGhvbGRvdXQgPT0gMSwgYmFnIDo9IGJhZ3MkYmFnc19uZXdkYXRhXQoKZGF0YV9yZWNvZGVkWywgLiguTiwgbWVhbihzdXJ2ZXllZCkpLCAuKGJhZyldW29yZGVyKGJhZyldCmBgYAoKIyMjIyBHZXQgbGFuZG1hcmsgcG9pbnRzCmBgYHtyfQpsYW5kbWFya3MgPSBnZXRMYW5kbWFya3MoZGF0YSA9IGRhdGFfcmVjb2RlZAogICAgICAgICAgICAgLCB2YXJzID0gdW5pcXVlKGNvdmFyc1tpbl9ib3RoID09IDEgfCBpbl9maWxlID09IDEsXSR2YXIpCiAgICAgICAgICAgICAsIG5fbGFuZG1hcmtzID0gMTIKICAgICAgICAgICAgICwgc3Vic2V0X2luZCA9IChkYXRhX3JlY29kZWQkdm90ZXJmaWxlID09IDEpKQoKWF9maWxlID0gbGFuZG1hcmtzJFhbZGF0YV9yZWNvZGVkJHZvdGVyZmlsZSA9PSAxLCBdClhfZmlsZV9ob2xkb3V0ID0gbGFuZG1hcmtzJFhbZGF0YV9yZWNvZGVkJGhvbGRvdXQgPT0gMSwgXQoKYGBgCgojIyMjIEVtYmVkIGZpbGUgaW4gZmVhdHVyZSBzcGFjZQoKQ2hvb3NlIHNjYWxlIHBhcmFtZXRlciAkXHNpZ21hJCBmb3IgUkJGIGtlcm5lbCAtIHVzZSBtZWRpYW4gaGV1cmlzdGljIGZvciBub3cKYGBge3J9CiMgcmJmMSA9IHJiZmRvdChzaWdtYSA9IDEpCiMgCiMgS19zaWdtYSA9IGtlcm5lbE1hdHJpeChyYmYxLCB4ID0gYXMubWF0cml4KFhfZmlsZSksIHkgPSBsYW5kbWFya3MkbGFuZG1hcmtzKQojIHNpZ21hID0gbWVkaWFuKEtfc2lnbWEpCgpgYGAKCkRvIGJhc2ljIERpc3QgUmVnCmBgYHtyfQp0cmFpbl9iYWcgPSBkYXRhX3JlY29kZWRbdm90ZXJmaWxlID09IDEsIGJhZ10KdGVzdF9iYWcgPSBkYXRhX3JlY29kZWRbaG9sZG91dCA9PSAxLCBiYWddCgojIGNhbGN1bGF0ZSBkZXBlbmRlbnQgdmFyIGluIGVhY2ggYmFnCllfc3Z5X2JhZyA9IGRhdGFfcmVjb2RlZFtzdXJ2ZXllZCA9PSAxLCAuKHlfbWVhbiA9IG1lYW4oeV9kZW0pKSwgYmFnXVtvcmRlcihiYWcpXQoKIyBnZXQgZmVhdHVyZXMKZmVhdHVyZXMgPSBnZXRGZWF0dXJlcyh0cmFpbiA9IFhfZmlsZQogICAgICAgICAgICAgICAgICAgICAgICwgdHJhaW5fYmFnID0gdHJhaW5fYmFnCiAgICAgICAgICAgICAgICAgICAgICAgLCB0ZXN0ID0gWF9maWxlX2hvbGRvdXQKICAgICAgICAgICAgICAgICAgICAgICAsIHRlc3RfYmFnID0gdGVzdF9iYWcKICAgICAgICAgICAgICAgICAgICAgICAsIGxhbmRtYXJrcyA9IGxhbmRtYXJrcyRsYW5kbWFya3MKICAgICAgICAgICAgICAgICAgICAgICAsIHNpZ21hID0gMC4wMSkKCiMgZG8gYmFzaWMgRFIKZml0X2Jhc2ljRFIgPSBmaXRMYXNzbyhtdV9oYXQgPSBmZWF0dXJlcyRtdV9oYXQKICAgICAgICAgICwgWV9iYWcgPSBZX3N2eV9iYWckeV9tZWFuCiAgICAgICAgICAsIHBoaV94X3RyYWluID0gZmVhdHVyZXMkcGhpX3hfdHJhaW4KICAgICAgICAgICwgcGhpX3hfdGVzdCA9IGZlYXR1cmVzJHBoaV94X3Rlc3QKICAgICAgICAgICkKCiMgc2NvcmUgdGhlIGZpbGUKZGF0YV9yZWNvZGVkW3ZvdGVyZmlsZSA9PSAxLCB5X2RlbV9iYXNpY2RyIDo9IGZpdF9iYXNpY0RSJFlfdHJhaW5dCmRhdGFfcmVjb2RlZFtob2xkb3V0ID09IDEsIHlfZGVtX2Jhc2ljZHIgOj0gZml0X2Jhc2ljRFIkWV90ZXN0XQoKZGF0YV9yZWNvZGVkWywgLk4sIC4oaXMubmEoeV9kZW1fYmFzaWNkcikpXQoKY2FsY01TRShZID0gZGF0YV9yZWNvZGVkJHlfZGVtLCBZX3ByZWQgPSBkYXRhX3JlY29kZWQkeV9kZW1fYmFzaWNkcikKYGBgCgoKYGBge3J9CgpkYXRhX3JlY29kZWRbdm90ZXJmaWxlID09IDEgfCBob2xkb3V0ID09IDEsIHlfZGVtX2Jhc2ljZHJfZGVjIDo9IGN1dCh5X2RlbV9iYXNpY2RyLCBicmVha3MgPSBxdWFudGlsZSh5X2RlbV9iYXNpY2RyLCBwcm9icyA9IHNlcSgwLDEsMC4xKSksIGxhYmVscyA9IDE6MTAsIGluY2x1ZGUubG93ZXN0ID0gVCldCgpnZ3Bsb3QoZGF0YV9yZWNvZGVkW3ZvdGVyZmlsZSA9PSAxLCAuKHBjdF95X2RlbSA9IG1lYW4oeV9kZW0pKSwgYnkgPSB5X2RlbV9iYXNpY2RyX2RlY10sIGFlcyh4ID0geV9kZW1fYmFzaWNkcl9kZWMsIHkgPSBwY3RfeV9kZW0pKSArIAogIGdlb21fYmFyKHN0YXQgPSAnaWRlbnRpdHknKSArCiAgZ2d0aXRsZSgiUGN0IERlbSBzdXBwb3J0ZXIgYnkgc2NvcmUgZGVjaWxlIC0gdm90ZXJmaWxlIikKZ2dwbG90KGRhdGFfcmVjb2RlZFtob2xkb3V0ID09IDEsIC4ocGN0X3lfZGVtID0gbWVhbih5X2RlbSkpLCBieSA9IHlfZGVtX2Jhc2ljZHJfZGVjXSwgYWVzKHggPSB5X2RlbV9iYXNpY2RyX2RlYywgeSA9IHBjdF95X2RlbSkpICsgCiAgZ2VvbV9iYXIoc3RhdCA9ICdpZGVudGl0eScpICsKICBnZ3RpdGxlKCJQY3QgRGVtIHN1cHBvcnRlciBieSBzY29yZSBkZWNpbGUgLSBob2xkb3V0IikKYGBgCgpIeXBlcnBhcmFtZXRlciB0dW5pbmcKYGBge3J9CmJhZ2dpbmdfdmFycyA9IGZpbGVfYW5kX3N1cnZleV92YXJzCnJlZ3Jlc3Npb25fdmFycyA9IHVuaXF1ZShjb3ZhcnNbaW5fYm90aCA9PSAxIHwgaW5fZmlsZSA9PSAxLF0kdmFyKQoKZml0X2Jhc2ljRFIgPSBkb0Jhc2ljRFIoZGF0YSA9IGRhdGFfcmVjb2RlZAogICAgICAgICAgICAgICAgICAgICAsIGJhZ2dpbmdfdmFycyA9IGZpbGVfYW5kX3N1cnZleV92YXJzCiAgICAgICAgICAgICAgICAgICAgICwgcmVncmVzc2lvbl92YXJzID0gdW5pcXVlKGNvdmFyc1tpbl9ib3RoID09IDEgfCBpbl9maWxlID09IDEsXSR2YXIpCiAgICAgICAgICAgICAgICAgICAgICwgb3V0Y29tZSA9ICd5X2RlbScKICAgICAgICAgICAgICAgICAgICAgLCBuX2JhZ3MgPSA3NQogICAgICAgICAgICAgICAgICAgICAsIG5fbGFuZG1hcmtzID0gMTAwCiAgICAgICAgICAgICAgICAgICAgICwgc2lnbWEgPSAwLjAxCiAgICAgICAgICAgICAgICAgICAgICwgYmFnZ2luZ19pbmQgPSAnc3VydmV5ZWQnCiAgICAgICAgICAgICAgICAgICAgICwgdHJhaW5faW5kID0gJ3ZvdGVyZmlsZScKICAgICAgICAgICAgICAgICAgICAgLCB0ZXN0X2luZCA9ICdob2xkb3V0JykKCgpuYW1lcyhmaXRfYmFzaWNEUikKYGBgCgoKYGBge3J9CnBhcmFtX3NpZ21hID0gZXhwKC1zZXEoMSw2LCBsZW5ndGgub3V0ID0gMTApKQpwYXJhbV9ubGFuZG1hcmtzID0gcm91bmQoZXhwKHNlcSgyLjUsIDcsIGxlbmd0aC5vdXQgPSAxMCkpKQpwYXJhbV9uYmFncyA9IGMoMjUsIDUwLCA3NSwgMTAwLCAxNTAsIDIwMCkKCnBhcmFtX2dyaWQgPSBleHBhbmQuZ3JpZChzaWdtYSA9IHBhcmFtX3NpZ21hLCBuX2JhZ3MgPSBwYXJhbV9uYmFncywgbl9sYW5kbWFya3MgPSBwYXJhbV9ubGFuZG1hcmtzKQoKZ3JpZF9zZWFyY2ggPSBhcHBseShwYXJhbV9ncmlkLCAxLCBmdW5jdGlvbihwKXsKICBjYXQocCkKICBjYXQoJ1xuJykKICBtc2UgPSB0cnlDYXRjaChkb0Jhc2ljRFIoZGF0YSA9IGRhdGFfcmVjb2RlZAogICAgICAgICAgICAgICAgICAgICAsIGJhZ2dpbmdfdmFycyA9IGZpbGVfYW5kX3N1cnZleV92YXJzCiAgICAgICAgICAgICAgICAgICAgICwgcmVncmVzc2lvbl92YXJzID0gdW5pcXVlKGNvdmFyc1tpbl9ib3RoID09IDEgfCBpbl9maWxlID09IDEsXSR2YXIpCiAgICAgICAgICAgICAgICAgICAgICwgb3V0Y29tZSA9ICd5X2RlbScKICAgICAgICAgICAgICAgICAgICAgLCBuX2JhZ3MgPSBwWzJdCiAgICAgICAgICAgICAgICAgICAgICwgbl9sYW5kbWFya3MgPSBwWzNdCiAgICAgICAgICAgICAgICAgICAgICwgc2lnbWEgPSBwWzFdCiAgICAgICAgICAgICAgICAgICAgICwgYmFnZ2luZ19pbmQgPSAnc3VydmV5ZWQnCiAgICAgICAgICAgICAgICAgICAgICwgdHJhaW5faW5kID0gJ3ZvdGVyZmlsZScKICAgICAgICAgICAgICAgICAgICAgLCB0ZXN0X2luZCA9ICdob2xkb3V0JwogICAgICAgICAgICAgICAgICApJG1zZV90ZXN0CiAgICAgICAgICAgICAgICAgLCBlcnJvciA9IGZ1bmN0aW9uKGUpIHJldHVybigwKSkKICAKICBtc2UKfSkKYGBgCgoKCmBgYHtyfQpnZ3Bsb3QoZGF0YS50YWJsZShwYXJhbV9ncmlkLCBncmlkX3NlYXJjaCksIGFlcyh4ID0gbl9sYW5kbWFya3MsIHkgPSBncmlkX3NlYXJjaCkpICsgZ2VvbV9wb2ludCgpCgpnZ3Bsb3QoZGF0YS50YWJsZShwYXJhbV9ncmlkLCBncmlkX3NlYXJjaCksIGFlcyh4ID0gbG9nKG5fbGFuZG1hcmtzKSwgeSA9IGxvZyhzaWdtYSkpKSArIGdlb21fdGlsZShhZXMoZmlsbCA9IC1sb2coZ3JpZF9zZWFyY2gpKSkgKyBmYWNldF9ncmlkKH5uX2JhZ3MpCmBgYAoK